{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# KKT Callbacks with IntPoint\n",
    "\n",
    "The interior point solver can sometimes be sped up dramatically by exploiting the simultanious structure of $Q$,$A$ and $G$. Here we will consider the simple problem:\n",
    "\n",
    "$$\\mbox{minimize}\\quad\\frac{1}{2}x^{T}Qx-b^{T}x,\\qquad\\mbox{s.t.}\\quad x\\geq0.$$\n",
    "\n",
    "The solver can be called with no special parameters,"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      " > INTERIOR POINT SOLVER v0.7 (July 2016)\n",
      "\n",
      "            Optimality                      Objective              Infeasibility       \n",
      "\n",
      "\u001b[1m   Iter  │  prFeas    duFeas    muFeas   │  pobj      dobj      │  icertp    icertd   │  refine   \u001b[0m\n",
      "      1  │  1.4e+00   4.3e+01   1.7e+00  │  -2.6e+01  -3.8e+03  │  NaN       1.0e+01  │  0\n",
      "      2  │  2.8e-01   8.4e+00   3.2e-01  │  -1.5e+01  -8.0e+02  │  NaN       2.0e+01  │  1\n",
      "      3  │  7.2e-02   2.2e+00   2.0e-01  │  -3.2e+00  -2.7e+02  │  NaN       6.8e+01  │  1\n",
      "      4  │  1.7e-02   5.3e-01   9.6e-02  │  -1.3e+00  -7.5e+01  │  NaN       1.5e+02  │  1\n",
      "      5  │  3.3e-03   1.0e-01   4.0e-02  │  -1.1e+00  -1.9e+01  │  NaN       2.5e+02  │  1\n",
      "      6  │  1.4e-14   5.9e-17   7.9e-03  │  -2.1e+00  -3.7e+00  │  NaN       2.9e+02  │  1\n",
      "      7  │  8.7e-15   3.4e-17   1.3e-03  │  -2.4e+00  -2.6e+00  │  NaN       2.5e+02  │  1\n",
      "      8  │  8.4e-15   3.1e-17   2.1e-04  │  -2.5e+00  -2.5e+00  │  NaN       2.4e+02  │  1\n",
      "      9  │  8.9e-15   3.4e-17   1.8e-05  │  -2.5e+00  -2.5e+00  │  NaN       2.4e+02  │  1\n",
      "     10  │  8.7e-15   3.5e-17   1.4e-06  │  -2.5e+00  -2.5e+00  │  NaN       2.4e+02  │  1\n",
      "     11  │  8.6e-15   3.6e-17   1.1e-07  │  -2.5e+00  -2.5e+00  │  NaN       2.4e+02  │  1\n",
      "\n",
      " > EXIT -- Below Tolerance!\n",
      "\n",
      "  5.597993 seconds (37.52 k allocations: 1.419 GB, 9.98% gc time)\n"
     ]
    }
   ],
   "source": [
    "using ConicIP\n",
    "\n",
    "n = 1000\n",
    "\n",
    "Q = sparse(randn(n,n)); Q = Q'*Q;\n",
    "c = ones(n,1);\n",
    "A = speye(n);\n",
    "b = zeros(n,1);\n",
    "𝐾 = [(\"R\",n)];\n",
    "\n",
    "@time conicIP( Q , c , A , b , 𝐾 , verbose = true);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The speed of the solver is reasonable, as the deault solver exploits the sparsity of the constraint matrix. We can do better, however.\n",
    "\n",
    "___\n",
    "\n",
    "## KKT Callbacks\n",
    "To speed up the solver, we require a function `kktsolver` which returns a function `solve3x3gen` which, in turn returns a function solving the KKT system\n",
    "$$\\left(\\begin{array}{ccc}\n",
    "Q & G^{T} & A^{T}\\\\\n",
    "G\\\\\n",
    "A &  & -F^{T}F\n",
    "\\end{array}\\right)\n",
    "\\left(\\begin{array}{c}\n",
    "a\\\\\n",
    "c\\\\\n",
    "b\n",
    "\\end{array}\\right) = \\left(\\begin{array}{c}\n",
    "x\\\\\n",
    "z\\\\\n",
    "y\n",
    "\\end{array}\\right)$$\n",
    "\n",
    "Where $F$ is a `Block` matrix with blocks corrosponding to the cones specified in 𝐾,\n",
    "$$F[i]=\\begin{cases}\n",
    "\\mbox{diag}(u_i) & \\mbox{if }K_{i}\\mbox{ is \"R\"}\\\\\n",
    "α_iJ+u_i u_i^T & \\mbox{if }K_{i}\\mbox{ is \"Q\"} \\quad\n",
    "\\mbox{ for some } u_i, \\alpha_i \\mbox{ $J$ is the hyperbolic identity }\\\\\n",
    "A_{U_i} & \\mbox{if }K_i\\mbox{ is \"S\"  }\\, \\quad \\mbox{where } A_{U_i}x=\\mbox{vec}(U_i^{T}\\mbox{mat}(x)U_i)  \\mbox{ for all $x$}.\n",
    "\\end{cases}\n",
    "$$\n",
    "\n",
    "The operation $vec$ is the vectorization operator on symmetric matrices \n",
    "$$\n",
    "\\mbox{vec}(U)=(U_{11},\\sqrt{2}U_{21},\\dots,\\sqrt{2}U_{p1},U_{22},\\sqrt{2}U_{32},\\dots,\\sqrt{2}U_{p2},\\dots,U_{p-1,p-1},\\sqrt{2}U_{p,p-1},U_{pp})\n",
    "$$\n",
    "\n",
    "and $\\mbox{mat}$ is the inverse transformation back to a symmetric matrix. The matrix $\\mbox{diag}(u)$ is represented by type `Diag`, $αI+uu$ is represented by type `SymWoodbury`, and  $A_{U_i}$ is represented by the type `VecCongurance.`\n",
    "\n",
    "In this example, since we have no linear constraints, $G$ is empty, and our KKT system is\n",
    "$$\n",
    "\\left(\\begin{array}{cc}\n",
    "Q & I\\\\\n",
    "I & -F^{T}F\n",
    "\\end{array}\\right)\\left(\\begin{array}{c}\n",
    "a\\\\\n",
    "b\n",
    "\\end{array}\\right)=\\left(\\begin{array}{c}\n",
    "x\\\\\n",
    "y\n",
    "\\end{array}\\right)\n",
    "$$\n",
    "\n",
    "The system can be solved by pivoting on the second block, as follows:\n",
    "$$\n",
    "\\big(Q+(F^TF)^{-1} \\big)\\,a=x+(F^TF)^{-1}y,\\qquad b=(F^TF)^{-1}(a-y)\n",
    "$$\n",
    "\n",
    "Because we only have polyhedral constraints, $F^{-2}$ is a diagonal matrix, thus the first equation is a diagonal perturbation to $Q$ which can be solved via a Cholesky Factorization. Pivoting allows us to solve a 1000x1000 system rather than a 2000x2000 system (albeit with a some sparsity structure)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      " > INTERIOR POINT SOLVER v0.7 (July 2016)\n",
      "\n",
      "            Optimality                      Objective              Infeasibility       \n",
      "\n",
      "\u001b[1m   Iter  │  prFeas    duFeas    muFeas   │  pobj      dobj      │  icertp    icertd   │  refine   \u001b[0m\n",
      "      1  │  1.3e+00   4.2e+01   1.8e+00  │  -2.2e+01  -3.6e+03  │  NaN       1.2e+01  │  0\n",
      "      2  │  2.4e-01   7.6e+00   2.8e-01  │  -1.3e+01  -6.9e+02  │  NaN       2.0e+01  │  1\n",
      "      3  │  7.3e-02   2.3e+00   1.8e-01  │  -6.0e+00  -2.5e+02  │  NaN       5.6e+01  │  1\n",
      "      4  │  2.0e-02   6.3e-01   9.7e-02  │  -2.4e+00  -8.3e+01  │  NaN       1.3e+02  │  1\n",
      "      5  │  3.0e-03   9.6e-02   3.8e-02  │  -1.3e+00  -1.8e+01  │  NaN       2.5e+02  │  1\n",
      "      6  │  1.2e-14   4.6e-17   7.5e-03  │  -2.1e+00  -3.6e+00  │  NaN       2.9e+02  │  1\n",
      "      7  │  8.8e-15   2.9e-17   1.0e-03  │  -2.4e+00  -2.6e+00  │  NaN       2.6e+02  │  1\n",
      "      8  │  8.2e-15   2.6e-17   1.4e-04  │  -2.4e+00  -2.5e+00  │  NaN       2.5e+02  │  1\n",
      "      9  │  8.5e-15   2.9e-17   1.3e-05  │  -2.4e+00  -2.4e+00  │  NaN       2.5e+02  │  1\n",
      "     10  │  8.9e-15   3.2e-17   9.8e-07  │  -2.4e+00  -2.4e+00  │  NaN       2.5e+02  │  1\n",
      "\n",
      " > EXIT -- Below Tolerance!\n",
      "\n",
      "  2.266916 seconds (37.09 k allocations: 1.035 GB, 13.42% gc time)\n"
     ]
    }
   ],
   "source": [
    "function kktsolver(Q, A, G, cone_dims)\n",
    "    \n",
    "    function solve3x3gen(F, F⁻¹)\n",
    "\n",
    "      invFᵀF = inv(F'F)\n",
    "      QpD⁻¹ = cholfact(Q + spdiagm( (F[1].diag).^(-2) ))\n",
    "\n",
    "      function solve3x3(x, z, y)\n",
    "\n",
    "        a = QpD⁻¹\\(x + A'*(invFᵀF*y))\n",
    "        b = invFᵀF*(y - A*a)\n",
    "        c = zeros(0,1)\n",
    "        return(a, c, b)\n",
    "\n",
    "      end\n",
    "\n",
    "    end\n",
    "    \n",
    "end\n",
    "\n",
    "@time sol = conicIP( Q , c , A , b , 𝐾 , kktsolver = kktsolver; verbose = true);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This results in a 5-fold improvement in speed, and a dramatic drop in memory usage!\n",
    "\n",
    "This pattern of pivoting on the third block happens often enough that we have encapsulated it in the convenience function `pivot`, which transforms a $2x2$ solver of the system\n",
    "\n",
    "$$\n",
    "\\left(\\begin{array}{cc}\n",
    "Q+A^{T}(F^{T}F)^{-1}A & G\\\\\n",
    "G & 0\n",
    "\\end{array}\\right)\\left(\\begin{array}{c}\n",
    "a\\\\\n",
    "b\n",
    "\\end{array}\\right) = \\left(\\begin{array}{c}\n",
    "x\\\\\n",
    "y\n",
    "\\end{array}\\right)\n",
    "$$\n",
    "\n",
    "into a $3x3$ solver. This is illustrated below"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      " > INTERIOR POINT SOLVER v0.7 (July 2016)\n",
      "\n",
      "            Optimality                      Objective              Infeasibility       \n",
      "\n",
      "\u001b[1m   Iter  │  prFeas    duFeas    muFeas   │  pobj      dobj      │  icertp    icertd   │  refine   \u001b[0m\n",
      "      1  │  1.3e+00   4.2e+01   1.8e+00  │  -2.2e+01  -3.6e+03  │  NaN       1.2e+01  │  0\n",
      "      2  │  2.4e-01   7.6e+00   2.8e-01  │  -1.3e+01  -6.9e+02  │  NaN       2.0e+01  │  1\n",
      "      3  │  7.3e-02   2.3e+00   1.8e-01  │  -6.0e+00  -2.5e+02  │  NaN       5.6e+01  │  1\n",
      "      4  │  2.0e-02   6.3e-01   9.7e-02  │  -2.4e+00  -8.3e+01  │  NaN       1.3e+02  │  1\n",
      "      5  │  3.0e-03   9.6e-02   3.8e-02  │  -1.3e+00  -1.8e+01  │  NaN       2.5e+02  │  1\n",
      "      6  │  1.3e-14   6.2e-17   7.5e-03  │  -2.1e+00  -3.6e+00  │  NaN       2.9e+02  │  1\n",
      "      7  │  9.2e-15   3.6e-17   1.0e-03  │  -2.4e+00  -2.6e+00  │  NaN       2.6e+02  │  1\n",
      "      8  │  9.2e-15   3.4e-17   1.4e-04  │  -2.4e+00  -2.5e+00  │  NaN       2.5e+02  │  1\n",
      "      9  │  8.8e-15   3.4e-17   1.3e-05  │  -2.4e+00  -2.4e+00  │  NaN       2.5e+02  │  1\n",
      "     10  │  8.7e-15   3.7e-17   9.8e-07  │  -2.4e+00  -2.4e+00  │  NaN       2.5e+02  │  1\n",
      "\n",
      " > EXIT -- Below Tolerance!\n",
      "\n",
      "  2.003137 seconds (48.64 k allocations: 1.036 GB, 9.69% gc time)\n"
     ]
    }
   ],
   "source": [
    "function kktsolver2x2(Q, A, G, cone_dims)\n",
    "    \n",
    "  function solve3x3gen(F, F⁻¹)\n",
    "\n",
    "    QpD⁻¹ = cholfact(Q + spdiagm( (F[1].diag).^(-2) ))\n",
    "    return (y, x) -> (QpD⁻¹\\y, zeros(0,1))\n",
    "\n",
    "  end\n",
    "    \n",
    "end\n",
    "\n",
    "@time sol = conicIP( Q , c , A , b , 𝐾 , kktsolver = pivot(kktsolver2x2); verbose = true);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And as a bonus we even get an extra boost in speed! \n",
    "___"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<style>\n",
       "\n",
       "/*div#maintoolbar, div#header {display: none !important;}\n",
       "*/\n",
       "\n",
       "#notebook {\n",
       "  padding-bottom: 20px\n",
       "}\n",
       "\n",
       "#notebook-container {\n",
       "font-family: Computer Modern;\n",
       "  width:900px;\n",
       "  padding: 20px;\n",
       "  background-color: #FFF;\n",
       "  min-height: 10;\n",
       "  border: 1px solid;\n",
       "}\n",
       "\n",
       "body {\n",
       "  background-color: white !important;\n",
       "}\n",
       "\n",
       ".MathJax_Display: {\n",
       "  margin: 100px\n",
       "}\n",
       "\n",
       "@font-face {\n",
       "  font-family: Computer Modern;\n",
       "src:url('http://gabgoh.github.io/cmunorm.ttf');\n",
       "}\n",
       "\n",
       "/* src: url('http://mirrors.ctan.org/\n",
       "fonts/cm-unicode/fonts/otf/cmunss.otf'); */\n",
       "\n",
       "div.cell{\n",
       "  width:850px;\n",
       "  margin-left:auto;\n",
       "  margin-right:auto;\n",
       "}\n",
       "\n",
       "div.output_subarea {\n",
       "    overflow-x: auto;\n",
       "    padding: 0.4em;\n",
       "    -webkit-box-flex: 1;\n",
       "    -moz-box-flex: 1;\n",
       "    box-flex: 1;\n",
       "    flex: 1;\n",
       "    max-width: calc(100%);\n",
       "}\n",
       "\n",
       ".output_wrapper {\n",
       "  padding-top : 10px;\n",
       "}\n",
       "\n",
       "div.input_area {\n",
       "    border: 0px solid #cfcfcf; \n",
       "    border-radius: 0px;\n",
       "    background: rgba(202, 202, 202, 0.28);\n",
       "    line-height: 1.21429em;\n",
       "    padding: 10px;\n",
       "}\n",
       "\n",
       "h1 { font-family: Computer Modern; }\n",
       "\n",
       ".prompt{ display:None; }\n",
       "\n",
       ".output_png img {\n",
       "    display: block !important;\n",
       "    margin-left: auto !important;\n",
       "    margin-right: auto !important;\n",
       "}\n",
       "\n",
       ".ui-wrapper {\n",
       "  margin-left:auto !important;\n",
       "  margin-right:auto !important\n",
       "}\n",
       "\n",
       "hr {\n",
       "    display: block;\n",
       "    height: 1px;\n",
       "    border: 0;\n",
       "    border-top: 1px solid #000;\n",
       "    margin: 1em 0;\n",
       "    padding: 0;\n",
       "}\n",
       "\n",
       ".err{\n",
       "  border : 0px !important;\n",
       "}\n",
       "\n",
       ".MathJax_Display{\n",
       "  padding-top: 12px;\n",
       "}\n",
       "\n",
       "\n",
       "</style>\n"
      ],
      "text/plain": [
       "HTML{ASCIIString}(\"<style>\\n\\n/*div#maintoolbar, div#header {display: none !important;}\\n*/\\n\\n#notebook {\\n  padding-bottom: 20px\\n}\\n\\n#notebook-container {\\nfont-family: Computer Modern;\\n  width:900px;\\n  padding: 20px;\\n  background-color: #FFF;\\n  min-height: 10;\\n  border: 1px solid;\\n}\\n\\nbody {\\n  background-color: white !important;\\n}\\n\\n.MathJax_Display: {\\n  margin: 100px\\n}\\n\\n@font-face {\\n  font-family: Computer Modern;\\nsrc:url('http://gabgoh.github.io/cmunorm.ttf');\\n}\\n\\n/* src: url('http://mirrors.ctan.org/\\nfonts/cm-unicode/fonts/otf/cmunss.otf'); */\\n\\ndiv.cell{\\n  width:850px;\\n  margin-left:auto;\\n  margin-right:auto;\\n}\\n\\ndiv.output_subarea {\\n    overflow-x: auto;\\n    padding: 0.4em;\\n    -webkit-box-flex: 1;\\n    -moz-box-flex: 1;\\n    box-flex: 1;\\n    flex: 1;\\n    max-width: calc(100%);\\n}\\n\\n.output_wrapper {\\n  padding-top : 10px;\\n}\\n\\ndiv.input_area {\\n    border: 0px solid #cfcfcf; \\n    border-radius: 0px;\\n    background: rgba(202, 202, 202, 0.28);\\n    line-height: 1.21429em;\\n    padding: 10px;\\n}\\n\\nh1 { font-family: Computer Modern; }\\n\\n.prompt{ display:None; }\\n\\n.output_png img {\\n    display: block !important;\\n    margin-left: auto !important;\\n    margin-right: auto !important;\\n}\\n\\n.ui-wrapper {\\n  margin-left:auto !important;\\n  margin-right:auto !important\\n}\\n\\nhr {\\n    display: block;\\n    height: 1px;\\n    border: 0;\\n    border-top: 1px solid #000;\\n    margin: 1em 0;\\n    padding: 0;\\n}\\n\\n.err{\\n  border : 0px !important;\\n}\\n\\n.MathJax_Display{\\n  padding-top: 12px;\\n}\\n\\n\\n</style>\\n\")"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "HTML(readall(open(\"style.css\")))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 0.4.5",
   "language": "julia",
   "name": "julia-0.4"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "0.4.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
